#include "base.h"

int main(int argc,char **argv)
{
          init_4F_projector();


	  //////////////////
	  // Input   ///////
	  //////////////////
          int Nr=2;
          char kappa_sea[1024];
          int conf_ini, conf_end, jump;
          int Nconf, Nex;
          char dir_input[1024];
          char dir_output[1024];

	  int beta100;
	
	  int imom;
	  int imass;
	  int irc; 

          char momentum_order[1024];
          char gamma[1024];
	

	  ////// Parsing command line

    	  for(int i=1; i<argc; i++) {
	    if(strncmp(argv[i],"-kappa_sea",10)==0)  sscanf(argv[i],"-kappa_sea=%s",kappa_sea);
	    if(strncmp(argv[i],"-conf_ini",9)==0)  sscanf(argv[i],"-conf_ini=%d",&conf_ini);
	    if(strncmp(argv[i],"-jump",5)==0)  sscanf(argv[i],"-jump=%d",&jump);
	    if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
	    if(strncmp(argv[i],"-Nex",4)==0)  sscanf(argv[i],"-Nex=%d",&Nex);
            if(strncmp(argv[i],"-dir_input",10)==0)  sscanf(argv[i],"-dir_input=%s",dir_input);
            if(strncmp(argv[i],"-dir_output",11)==0)  sscanf(argv[i],"-dir_output=%s",dir_output);
	    if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T); 
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L); 
            if(strncmp(argv[i],"-beta100",8)==0)  sscanf(argv[i],"-beta100=%d",&beta100);

	    if(strncmp(argv[i],"-imass",6)==0)  sscanf(argv[i],"-imass=%d",&imass); 
            if(strncmp(argv[i],"-imom",5)==0)  sscanf(argv[i],"-imom=%d",&imom);
	    if(strncmp(argv[i],"-irc",4)==0) sscanf(argv[i],"-irc=%d",&irc);

             if(strncmp(argv[i],"-momentum_order",15)==0)  sscanf(argv[i],"-momentum_order=%s",momentum_order);
             if(strncmp(argv[i],"-gamma",6)==0)  sscanf(argv[i],"-gamma=%s",gamma);

    	  }

          init_base_gamma(gamma);

        ///////////////////////////////
        //// r-combinations
        ///////////////////////////////

	// Como:
	// OpmPX(0010)=OpmPX(1011)=OpmPX(1110)*=OpmPX(1000)*
	// OpmPX(0001)=OpmPX(1101)=OpmPX(0111)*=OpmPX(0100)*
	// (Comprobado numéricamente)
	// (Parece lógico porque es cambiar la parte TW - OS pero buscar la simetria) //TODO
	// Basta calcular
	// OpmPX(0001) y OpmPX(1110) y al promediar quedarnos con la parte real

          int *rvector=malloc(sizeof(int)*(4)); whichr_sort(irc, rvector);
          int r1,r2,r3,r4; r1=rvector[0]; r2=rvector[1]; r3=rvector[2]; r4=rvector[3];

	  //////////////////////////////
	  // print some input
	  ///////////////////////////////
	  Nconf=Nclust*Nex;
	  conf_end=conf_ini+jump*(Nconf-1);
          printf("=================================================================================================================\n");
	  printf("------> imom=%d <-------------\n",imom);
	  printf("------> imass=%d <-------------\n",imass);
	  printf("------> r=(%d,%d,%d,%d) <---------\n",r1,r2,r3,r4);
          printf("conf_ini=%d, conf_end=%d jump=%d\nNconf=%d, Nclust=%d, Nex=%d\n",conf_ini, conf_end,jump,Nconf,Nclust,Nex);
          printf("imass=%d\n",imass);
          printf("L=%d T=%d\n",L,T);
          printf("Nr=%d \n",Nr);
          printf("dir_input=%s \n",dir_input);
          printf("dir_output=%s \n",dir_output);
          printf("===================================================================================================================\n");
	


   /////////////////////////
  //// Momenta set
  //////////////////////////////
  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom_%d.dat",beta100);
  read_momenta_data(momentafile);
  double *nn=(double*)malloc(sizeof(double)*(Nmom));
  compute_nn(momentum_order,nn);
  Momentum(momentum_order);

   ///////////////////
   // IO files
   //////////////////

          char s0ff_base_input[1024];
          char s0ft_input[1024];
          char s0ff_base_output[1024];
          char z_base_output[1024];
          char Lambda_base_output[1024];
	  char file[1024];
 
          sprintf(s0ff_base_input,"%s/s0ft-mom-m/s0ft",dir_input);
          sprintf(z_base_output,"%s/Z_ALLP/Z",dir_output);
          sprintf(s0ff_base_output,"%s/FF_ALLP/ff",dir_output);
          sprintf(Lambda_base_output,"%s/LAMBDA_ALLP/Lambda",dir_output);


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

           /////////////////////////////
	   // Read Fourier transform
	   //////////////////////////////
	
          //def
	  colorcolorspinspin **s0FT_r=malloc(sizeof(colorcolorspinspin*)*(Nconf));
	  for (int iconf=0; iconf<Nconf; iconf++)   s0FT_r[iconf]=malloc(sizeof(colorcolorspinspin)*(Nr));
	
	  int conf_to_read;
	  FILE *fin=NULL;


           printf("====================\n");
           printf("READING...\n");

           for (int r=0; r<Nr; r++)
	   {
           for (int iconf=0; iconf<Nconf; iconf++)
 	   {
	           conf_to_read=conf_ini+iconf*jump;
                   //printf("  ... conf %d\n",conf_to_read);
                   sprintf(s0ft_input,"%s%d_%s_%04d_%04d_m%d.out",s0ff_base_input,r,kappa_sea,conf_to_read,imom,imass);
                   fin=fopen(s0ft_input,"rb+");
                   if(fin == (FILE *)NULL) {
                   sprintf(s0ft_input,"%s%d__%04d_%04d_m%d.out",s0ff_base_input,r,conf_to_read,imom,imass);
                   fin=fopen(s0ft_input,"rb+");
			if(fin == (FILE *)NULL) {
                        fprintf(stderr,"Unable to open file %s for reading\n",s0ft_input);
                        return EXIT_FAILURE;
			}
                   }
                  // printf("im=%d\n",imass);
			for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
			fread(&s0FT_r[iconf][r][ic_int][ic_ext][is_int][is_ext],sizeof(complex),1,fin);
		   fclose(fin);
    	
   	   }
	   }
	   printf("====================\n\n\n");


	    ////////////////////////////////////
            /////         CLUSTERIZE      //////
	    ////////////////////////////////////
	    printf("====================\n");
            printf("CLUSTERIZE\n");

	   complex   aux_conf[Nconf], aux_clust[Nclust+1];

          //def
          colorcolorspinspin **s0FT=malloc(sizeof(colorcolorspinspin*)*(Nclust+1));
	  for (int iclust=0; iclust<Nclust+1; iclust++)  s0FT[iclust]=malloc(sizeof(colorcolorspinspin)*(Nr));

	  for (int r=0; r<Nr; r++)
		for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
		{	
	               for(int iconf=0;iconf<Nconf;iconf++) aux_conf[iconf]=s0FT_r[iconf][r][ic_int][ic_ext][is_int][is_ext];
                       conf2clust(aux_conf, aux_clust,Nconf,Nclust,Nex);
                       for(int iclust=0;iclust<Nclust+1;iclust++) s0FT[iclust][r][ic_int][ic_ext][is_int][is_ext]=aux_clust[iclust];
		}
	   
            printf("====================\n\n\n");




//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

	/////////////////////////////////////////////////////
        // FORM FACTOR OF THE FOURIER-TRANSFORMED PROPAGATOR
	//////////////////////////////////////////////////////

        printf("====================\n");
        printf("FORM FACTOR\n");

	//def
        complex ***s0_ff=malloc(sizeof(complex**)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++)  s0_ff[iclust]=malloc(sizeof(complex*)*(Nr));
        for (int iclust=0; iclust<Nclust+1; iclust++)  for (int r=0; r<Nr; r++) s0_ff[iclust][r]=malloc(sizeof(complex)*(3));
	
		        for (int iclust=0; iclust<Nclust+1; iclust++) for (int r=0; r<Nr; r++)
			compute_s0_ff(s0FT[iclust][r],s0_ff[iclust][r],SinP[imom],SinP2[imom]);

      

        printf("====================\n\n\n");
        

          
	////////////////////////////////
	// INVERTING		////////
	////////////////////////////////

        printf("====================\n");
	printf("INVERTING PROPAGATOR\n");

        colorcolorspinspin **s0FT_inv=malloc(sizeof(colorcolorspinspin*)*(Nclust+1));
	for (int iclust=0; iclust<Nclust+1; iclust++)  s0FT_inv[iclust]=malloc(sizeof(colorcolorspinspin)*(Nr));



	complex **mdum=(complex**)malloc(sizeof(complex)*(12));
	for(int i=0; i<12; i++) mdum[i]=(complex*)malloc(sizeof(complex)*(12));

        complex **idum=(complex**)malloc(sizeof(complex)*(12));
        for(int i=0; i<12; i++) idum[i]=(complex*)malloc(sizeof(complex)*(12));

	colorcolorspinspin s0_dump;

	int ind_ext, ind_int;
	int iflag=0; 

	for (int r=0; r<Nr; r++)
        
        {
	    for (int iclust=0; iclust<Nclust+1; iclust++) 
            {
                        for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
                         {
				ind_ext=is_ext+4*(ic_ext);
				ind_int=is_int+4*(ic_int);
				mdum[ind_int][ind_ext]=s0FT[iclust][r][ic_int][ic_ext][is_int][is_ext];			

			}
             		inverse(mdum,12,idum,iflag);
			if (iflag==1) {break; printf("Light quark propagator is singular\n");}
			for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
                         {
                                ind_ext=is_ext+4*(ic_ext);
                                ind_int=is_int+4*(ic_int);
				s0FT_inv[iclust][r][ic_int][ic_ext][is_int][is_ext]=idum[ind_int][ind_ext];
			}
			/////////////////////////
			// Inversion Test
			///////////////////////
			for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
			{
				s0_dump[ic_int][ic_ext][is_int][is_ext]=0.;
				for (int ic_dump=0; ic_dump<3; ic_dump++) for (int is_dump=0; is_dump<4; is_dump++)
						s0_dump[ic_int][ic_ext][is_int][is_ext]=s0_dump[ic_int][ic_ext][is_int][is_ext]+s0FT[iclust][r][ic_int][ic_dump][is_int][is_dump]*s0FT_inv[iclust][r][ic_dump][ic_ext][is_dump][is_ext];
			}
			for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
			{
				if (ic_ext != ic_int || is_ext != is_int) 
				{
					if (cabs(s0_dump[ic_int][ic_ext][is_int][is_ext])> 0.1) 
						printf("Inversion failed for imom=%d (s0*s0inv)_%d%d%d%d=%8.5f+I%8.5f \n",imom,ic_int,ic_ext,is_int,is_ext,creal(s0_dump[ic_int][ic_ext][is_int][is_ext]),cimag(s0_dump[ic_int][ic_ext][is_int][is_ext]));
				}
				else
                                 {
				       if (cabs(s0_dump[ic_int][ic_ext][is_int][is_ext])-1> 0.1) 
						printf("Inversion failed  for imom=%d (s0*s0inv)_%d%d%d%d=%8.5f+I%8.5f \n",imom,ic_int,ic_ext,is_int,is_ext,creal(s0_dump[ic_int][ic_ext][is_int][is_ext]),cimag(s0_dump[ic_int][ic_ext][is_int][is_ext]));
				}
				
			}
                        ////////////////////////////////////
                        //  end inversion test
                        ////////////////////////////////////

	     }
        }
	printf("====================\n\n\n");

       /////////////////////////////////////////////////////
        // FORM FACTOR OF THE INVERTED FOURIER-TRANSFORMED PROPAGATOR
        //////////////////////////////////////////////////////

	printf("====================\n");
        printf("FORM FACTOR OF THE INVERTED PROPAGATOR\n");
        complex ***s0_inv_ff=malloc(sizeof(complex**)*(Nclust+1));
	for (int iclust=0; iclust<Nclust+1; iclust++)  s0_inv_ff[iclust]=malloc(sizeof(complex*)*(Nr));
        for (int iclust=0; iclust<Nclust+1; iclust++)  for (int r=0; r<Nr; r++) s0_inv_ff[iclust][r]=malloc(sizeof(complex)*(3));


       for (int r=0; r<Nr; r++)
       
       {
		for (int iclust=0; iclust<Nclust+1; iclust++) 
                compute_s0_inv_ff(s0FT_inv[iclust][r],s0_inv_ff[iclust][r],SinP[imom],SinP2[imom],nn[imom]);
       }



    


	printf("====================\n\n\n");
  
        /////////////////////////////////////////////////////
        // Z_q=-i/12 * Tr(sum_mu gamma_mu*SinP(mu)*S^-1(p))/SinP2
        //////////////////////////////////////////////////////

        printf("====================\n");
        printf("Z_q\n");


        complex **Zq=malloc(sizeof(complex*)*(Nclust+1));
	for (int iclust=0; iclust<Nclust+1; iclust++)  Zq[iclust]=malloc(sizeof(complex)*(Nr));


        for (int r=0; r<Nr; r++)
        for (int iclust=0; iclust<Nclust+1; iclust++) 
        {
                        Zq[iclust][r]=s0_inv_ff[iclust][r][0];
        }






	//////////////////////////////////////////////////
        ////////////////////////////////////////////////////
       printf("====================\n\n\n");
       printf(" WRITING FF & ZQ OUTPUT IN %s\n",dir_output);

          
      for (int sigma=0; sigma<3; sigma++)
      {

        for (int r=0; r<Nr; r++)
        {
        // ff S(p)

          for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=s0_ff[iclust][r][sigma];
         sprintf(file,"ff_r%d_m%d_sigma%d",r,imass,sigma+1);
         write_output_onemomenta(dir_output,file,aux_clust,SinP2,imom);


        }
      }
      for (int sigma=0; sigma<3; sigma++)
      {

        for (int r=0; r<Nr; r++)
        {

        // ff S^-1(p)
         for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=s0_inv_ff[iclust][r][sigma];
         sprintf(file,"ff_inv_r%d_m%d_sigma%d",r,imass,sigma+1);
         write_output_onemomenta(dir_output,file,aux_clust,SinP2,imom);

        }
      }

    // Zq
    for (int r=0; r<Nr; r++)
     {
         for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=Zq[iclust][r];
         sprintf(file,"Zq_r%d_m%d",r,imass);
         write_output_onemomenta(dir_output,file,aux_clust,SinP2,imom);
     }


	free(Zq);
	free(s0_inv_ff); free(s0_ff);

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

	/////////////////////////////////////////////////////
	/////////////////////////////////////////////////////
	///////////////  Z4F  ///////////////////////////////
	//////////////////////////////////////////////////////
	//////////////////////////////////////////////////////

        printf("====================\n");
        printf(" Z4F\n");

        gamma_ccss *GOS=malloc(sizeof(gamma_ccss)*(Nconf));
        gamma_ccss *GTW=malloc(sizeof(gamma_ccss)*(Nconf));
      
	colorcolorspinspin S_xy[2];
	colorcolorspinspin S_yx[2];
	colorcolorspinspin Sinv_xy[2];
	colorcolorspinspin Sinv_yx[2];

	O5_ccccssss *OpPC_r=malloc(sizeof(O5_ccccssss)*(Nconf));
        O5_ccccssss *OpPC=malloc(sizeof(O5_ccccssss)*(Nclust+1));


	// Dynamical matrix
        O55  *DpPC=malloc(sizeof(O55)*(Nclust+1));


	printf("Computing non Amputed Green Function...\n");
        for (int iconf=0; iconf<Nconf; iconf++)
        {
	//	if(Nconf>80)
	//	{
	//	for (int part=0; part<80; part++) if (iconf==part*Nconf/80) printf("\n");
        //      printf("%d ... ",iconf);
	//	}
	
        ////////////////////////////////////////////////////////
        // Compute 2F NON amputed Green Function for each 
        // Gamma=I, g1, g2, g3, g0, g5, sigma23, sigma31, sigma12, 
        //        sigma01, sigma02, sigma03, g1g5, g2g5, g3g5, g0g5
        // G=\psi_1 Gamma \psi_2
        // with r1r2=00,01,10,11
        //////////////////////////////////////////////////////

         for (int igamma=0; igamma<16; igamma++)
         {
			// 1. Initialize
			for (int r=0; r<2; r++) 
			for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++) 
                        for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++) 
			{
				S_xy[r][ic1][ic2][is1][is2]=s0FT_r[iconf][r][ic1][ic2][is1][is2];
				S_yx[r][ic1][ic2][is1][is2]=((base_gamma[5]).entr[is1])*conj(s0FT_r[iconf][!r][ic2][ic1][(base_gamma[5]).pos[is2]][(base_gamma[5]).pos[is1]])*(base_gamma[5]).entr[is2]; // Note the change in r 
			}
			// 2. Compute
                       for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++)
                       for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++)
			{
				GOS[iconf][igamma][ic1][ic2][is1][is2]=0.;
				GTW[iconf][igamma][ic1][ic2][is1][is2]=0.;
				for (int ic3=0; ic3<3 ; ic3++) for (int is3=0; is3<4 ; is3++){
					GOS[iconf][igamma][ic1][ic2][is1][is2]=GOS[iconf][igamma][ic1][ic2][is1][is2]+S_xy[r1][ic1][ic3][is1][is3]*((base_gamma[igamma]).entr[is3])*S_yx[r2][ic3][ic2][(base_gamma[igamma]).pos[is3]][is2];
                                        GTW[iconf][igamma][ic1][ic2][is1][is2]=GTW[iconf][igamma][ic1][ic2][is1][is2]+S_xy[r3][ic1][ic3][is1][is3]*((base_gamma[igamma]).entr[is3])*S_yx[r4][ic3][ic2][(base_gamma[igamma]).pos[is3]][is2];
				}

			}
       } // end igamma

     /////////////////////////////////////////////////////
     // Compute 4F NON Amputed Green Function from 2F
     // OpmPC_r OpmPV_r:
     // OpmPC_r=VV+AA,VV-AA,SS-PP,SS+PP,TT
     // OpmPV_r=VA+AV,VA-AV,SP-PS,SP+PS,TT~
     // Opm corresponds to eq.25 in bibbia while operators are defined in eq.12 and eq.16 coming from eq.7
     // Non Amputed Green Function are computed for each r1=r2=r3=-r4 type combination
     /////////////////////////////////////////////////////

        compute_OpPC( OpPC_r[iconf],  GOS[iconf], GTW[iconf]);


   } // end iconf


     /////////////////////
     // CLUSTERIZE
     ///////////////////
     // On one single configuration the quark propagator is not transationally invariant. It is only recovered after averaging over the configurations
     // Then we compute G over the initial configurations 
     // and then we average
     // G(p)=\sum_n G_n(p)/N
       printf("Clusterize...\n");
      for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++) for (int ic3=0; ic3<3 ; ic3++) for (int ic4=0; ic4<3 ; ic4++)
      for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++) for (int is3=0; is3<4 ; is3++) for (int is4=0; is4<4 ; is4++)
      for (int iop=0; iop<5; iop++)
      {
                       for(int iconf=0;iconf<Nconf;iconf++) aux_conf[iconf]=OpPC_r[iconf][iop][ic1][ic2][ic3][ic4][is1][is2][is3][is4];
                       conf2clust(aux_conf, aux_clust,Nconf,Nclust,Nex);
                       for(int iclust=0;iclust<Nclust+1;iclust++) OpPC[iclust][iop][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=aux_clust[iclust];
      }

    ////////////////////////////////////
    // Dynamical matrix:
    // D=Trace [ Proj * Lambda ]
    // Lambda= S^-1 S^-1 G S^-1 S^-1
    /////////////////////////////////////////

    printf("Computing Dynamical matrix\n");
    for (int iclust=0; iclust<Nclust+1; iclust++)
    {
                printf("%d ...\n",iclust);

       // Inverted propagator
        for (int r=0; r<2; r++)
         for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++)
          for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++)
          {
           Sinv_xy[r][ic1][ic2][is1][is2]=s0FT_inv[iclust][r][ic1][ic2][is1][is2];
           Sinv_yx[r][ic1][ic2][is1][is2]=((base_gamma[5]).entr[is1])*conj(s0FT_inv[iclust][!r][ic2][ic1][(base_gamma[5]).pos[is2]][(base_gamma[5]).pos[is1]])*(base_gamma[5]).entr[is2];
          }

	// Tree level projection
	// SS PP VV AA TT
	// SP PS VA AV TT5
	
	Projection_TL_OpPC(OpPC[iclust],Sinv_xy,Sinv_yx,r1,r2,r3,r4);	

	// Dynamical matrix

             for (int iop=0; iop<5; iop++)
             {

                DpPC[iclust][iop][0]=PpPC[0][0]*(VV[iop]+AA[iop]);
                DpPC[iclust][iop][1]=PpPC[1][1]*(VV[iop]-AA[iop]) + PpPC[1][2]*(SS[iop]-PP[iop]);
                DpPC[iclust][iop][2]=PpPC[2][1]*(VV[iop]-AA[iop]) + PpPC[2][2]*(SS[iop]-PP[iop]);
                DpPC[iclust][iop][3]=PpPC[3][3]*(SS[iop]+PP[iop]) + PpPC[3][4]*TT[iop];
                DpPC[iclust][iop][4]=PpPC[4][3]*(SS[iop]+PP[iop]) + PpPC[4][4]*TT[iop];
                
             }

    } // end iclust




    printf("====================\n\n\n");





///////////////////

      printf(" WRITING OUTPUT IN %s\n",dir_output);



   //////////////////////
   // DpmPC
   //////////////////


     // write .dat file with all D_ij and its error
     sprintf(file,"%s/%d%d%d%d/p%03d/DpPC_m%d.dat",dir_output,r1,r2,r3,r4,imom,imass);
     FILE *fout=fopen(file,"w");
     if(fout == (FILE *)NULL) {
              fprintf(stderr,"Unable to open file %s for writing\n",file);
              return EXIT_FAILURE;
     }
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
      {
                 for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]= DpPC[iclust][iop1][iop2];
                 fprintf(fout,"%8.10f %8.10f %8.10f  \n",SinP2[imom],creal(DpPC[Nclust][iop1][iop2]),creal(error_jack(aux_clust,Nclust)));
      }
      fclose(fout);

      // write each D_ij in different file for all cluster
      char dir_output_r[1024];
      sprintf(dir_output_r,"%s/%d%d%d%d/",dir_output,r1,r2,r3,r4);
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
      {
		sprintf(file,"DpPC%d%d_m%d",iop1+1,iop2+1,imass);
		for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=DpPC[iclust][iop1][iop2];
		write_output_onemomenta(dir_output_r,file,aux_clust,SinP2,imom);
      }

	printf("Calc_Z4F.c finished OK ! \n");

	free(GTW); free(GOS);  free(OpPC); free(OpPC_r) ; free(DpPC);
	
  return 0;
}
